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Abstract 

We study the problem of critical slowing-down for gauge-fixing algorithms (Landau 
gauge) in SU{2) lattice gauge theory on a 2-dimensional lattice. We consider five such 
algorithms, and lattice sizes ranging from 8^ to 36^ (up to 64? in the case of Fourier 
acceleration). We measure four different observables and we find that for each given 
algorithm they all have the same relaxation time within error bars. We obtain that: the 
so-called Los Alamos method has dynamic critical exponent z ~ 2, the overrelaxation 
method and the stochastic overrelaxation method have z ~ 1, the so-called Cornell 
method has z slightly smaller than 1 and the Fourier acceleration method completely 
eliminates critical slowing-down. A detailed discussion and analysis of the tuning of 
these algorithms is also presented. 



1 Introduction 

The lattice formulation of QCD provides a regularization which makes the gauge group 
compact, so that the Gibbs average of any gauge-invariant quantity is well-defined and thus 
gauge fixing is, in principle, not needed. However, to better understand the relationship 
between continuum and lattice models, one is led to consider gauge- dependent quantities on 
the lattice as well, which requires gauge fixing. It is therefore important to devise numerical 
algorithms to efficiently gauge fix a lattice configuration. The efficiency of these algorithms 
is even more important if the problem of existence of Gribov copies in the lattice is taken 
into account In fact, since usually it is not clear how an algorithm selects among 

different Gribov copies, numerical results using gauge fixing could dependQ on the gauge- 
fixing algorithm, making their interpretation conceptually difficult. In these cases, in order 
to analyze the dependence on the Gribov ambiguity p, ^, ^, several Gribov copies of the 
same thermalized configuration have to be produced, and therefore gauge fixing is extensively 
used. 



In this paper we study the problem of fixing the standard lattice Landau gauge condition 
P, 1^. As we will see in the next section, this condition is formulated as a minimization 
problem for the energy of a nonlinear cr-model with disordered couplings. Two basic deter- 
ministic local algorithms have been introduced to achieve this goal. Following we will 
refer0 to them as the "Los Alamos" method [|l^, |13| and the "Cornell" method [|l4l. Both 



methods are expected to perform poorly 0, especially as the volume of the lattice increases, 
due to the phenomenon of critical slowing-down, which afflicts Monte Carlo simulations of 
critical phenomena as well as some deterministic iterative methods, such as our gauge fix- 
ing.0 Roughly speaking, the problem is that, since the updating is local, the "information" 
travels at each step only from one lattice site to its nearest neighbors, executing a kind of 
random walk through the lattice; as a result, in order to get any significant change in con- 
figuration, we must wait a time of the order of the square of the typical "physical length" 
of the system, which is in our case the lattice size A^. More precisely, the relaxation time 
(measured in sweeps) for conventional local algorithms diverges as the square of the linear 
size of the system, or equivalently these methods have dynamic critical exponent z = 2. This 
A^^ behavior is not a serious difficulty for small lattices, and other aspects of the algorithms 
may be of greater importance in this case; but as one deals with progressively larger lattices 

^ Of course, this does not apply to gauge-invariant quantities. 

^ Some of the algorithms we consider here are also well-known for other numerical problems, and are 
usually referred to with other names. In particular, if we consider the problem of solving a linear system of 
equations |^ (or the equivalent problem of minimizing the related quadratic action), then the Los Alamos 
method corresponds to a non-linear version of the Gauss-Seidel method, while the overrelaxation method 
- and the Cornell method (see Sections 5 and 7) — correspond to non-linear versions of the successive 
overrelaxation method. On the other hand, from the point of view of minimizing a function, the Cornell 
method is a steepest- descent method [lO| , since the function is minimized in the direction of the local downhill 
gradient, while the Los Alamos method brings the "local" function to its unique absolute minimum. In both 
cases, the idea is to decrease the value of the minimizing function monotonically, namely these are descent 
methods |jll|. 

^ For an excellent introduction to the problem of critical slowing-down in Monte Carlo simulations, and 
also to some deterministic examples, see [[L5[. 
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(in order to approach the continuum hmit) this factor constitutes a severe hmitation. 

To overcome this problem, two solutions are available: one can either modify the local 
update, in such a way that the "length" of the move in the configuration space is increased 
|T^, and therefore this space is explored faster, or introduce some kind of global updating. 



in order to speed up the relaxation of the long-wavelength modes, which are the slow ones. 
Various methods, based on these two ideas, have been proposed: overrelaxation , stochas- 



tic overrelaxation ||T3|, multigrid schemes ||T3|, |T^, [T^, [T^, Fourier acceleration [|T4|, wavelet 



acceleration [^], etc. By analogy with other deterministic problems |TT| and with 

Monte Carlo simulations the "improved" local algorithms (such as overrelaxation 

and stochastic overrelaxation) are expected to reduce but not eliminate critical slowing- 
down (2; 1, as opposed to z ^ 2 for conventional methods), while global methods (such 
as multigrid, Fourier acceleration, etc.) hope to eliminate critical slowing-down completely 
{i.e., z = 0). In any case, a precise determination of the dynamic critical exponents for the 
different algorithms is of great importance, as are analyses and comparisons between the 
methods applied to the problem at hand. Such analyses have been partially done in some 
of the references given here, but we feel that systematic studies of these methods for the 
specific problem of Landau gauge fixing are lacking, especially with regard to the evaluation 
of dynamic critical exponents, although some of the algorithms we consider were extensively 
analyzed when applied to other numerical problems. The only algorithm thoroughly studied 
in the past was the multigrid [0, |13|, |18|, |19| and we will not consider it here. For our research. 



we have decided to study, besides the two "basic" algorithms (Los Alamos and Cornell), the 
standard overrelaxation and the stochastic overrelaxation (both applied to the Los Alamos 
method), which are very appealing for their simplicity and almost absence of overhead. We 
will also study the Fourier acceleration (applied to the Cornell method), for which not so 
many studies have been done up to now, and which is claimed to be the best method 
available today for Landau gauge fixing. 
In this work our goals are: 

1. studying the critical slowing-down for the various algorithms and finding accurate 
values for their dynamic critical exponents, 

2. analyzing the relative size of several quantities, usually employed in the literature to 
test the convergence of the gauge fixing, in order to understand which of them should 
be used in practical computations, 

3. doing, when necessary, a careful tuning of the algorithms, checking the "empirical" 
formulae commonly used for the optimal choice of the parameters, or finding a simple 
prescription when this formula is not known, 

4. comparing the computational costs of the algorithms, 

5. doing a simple analysis of the algorithms in order to get at least an idea of how they 
deal with the problem of critical slowing-down. 

Regarding the last point, we will essentially try to review what is known about the algorithms 
under analysis. Here, in fact, we are unable to do a more rigorous analysis (in the style of 
2^ or ||2^ for the gaussian model), due to the non-linearity of the update and the presence 



of "random" link- dependent coupling constants. 
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The paper is organized as follows. In Section 2 we present a pedagogical review of the 
problem of Landau gauge fixing on the lattice. In particular, to make the analysis of the 
efficiency of the local algorithms a little more quantitative, we introduce two functions: the 
variation of the minimizing function at each site A£{x) , and the "length" of the update 
"Dix) , interpreted as a move in the configuration space {g{x)} . In Section 3 we define the 
update for the various algorithms and we explain how each one fights the problem of critical 
slowing-down. The quantities for which the relaxation time r will be measured are introduced 
in Section 4. In Section 5 the problem of the tuning of the algorithms is addressed, and we 
try a quantitative analysis to find simple formulae for the optimal choice of the parameters. 
Finally, in Section 6, we give some more details about the numerical simulations and the 
computational aspects of our work and, in Section 7, we present our results and conclusions. 

The main difficulty we had to overcome in this project was the severe lack of computer 
time, which restricted us to dealing with the SU (2) gauge group, instead of the more inter- 
esting SU{?)) case, and with small lattice sizes. On the other hand, the use of SU (2) makes 
the analysis of the algorithms simpler, and with the values of (3 that we consider (see Section 
4) no significant finite-size corrections are expected to occur and the use of small lattices is 
justified. 

A further difficulty is the definition of constant physics, necessary for finding the dynamic 
critical exponents that characterize each algorithm (see Section 4). This definition is very 
simple only in dimension d = 2 and this is the case we will consider here, leaving the extension 
of this work to four dimensions to a future paper . 



Nevertheless, we believe that this work presents a comprehensive analysis and comparison 
of the different methods considered, and enough evidence for the evaluation of their dynamic 
critical exponents. Our findings for the exponents are basically confirmations of what is 
generally accepted, with the exception of the value slightly smaller than one for the Cornell 
method, a fact that we try to interpret in Sections 5 and 7. 

The total computer time used in our simulations was about 225 hours on an IBM RS- 
6000/250 machine. 

2 Landau Gauge-Fixing Condition on the Lattice 

Let us consider a standard Wilson action for SU{2) lattice gauge theory in d dimensions 



m ■■ 

S{{U})^^-Y. Y.\^-l^Tr[u,{x)U,{x + ae,)U-^\x + ae,)U;;\x)\ \ (2.1) 

where Ufj,{x) G SU{2) are the link variables, go is the bare coupling constant, a is the lattice 
spacing and is a unit vector in the positive fi direction. Sites are labeled by (i-dimensional 
vectors x . The lattice size in the /z direction is = a A*"^ where A*"^ is an integer. We 
assume periodic boundary conditions, i.e. x + L^e^ = x, and the lattice volume is given by 

d d 

V=l[L^ = a''l[N^. (2.2) 

fj,=i 11=1 
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The gauge field is defined as 



2ago 



(2.3) 



this variable approaches the classical gauge potential in the continuum limit. 
To fix the Landau gauge we look for a local minimum^ of the function 



£{{9}) 



l-^EE^Tr[f/(^)(a.) + [/(^)t(a.) 



(2.4) 



keeping the configuration {f7^(a;)} fixed. Here each g{x) G SU{2) is a site variable, Q = 
{g(x)} is a gauge transformation, and U^^\x) is given by 



U,{x) ^ Ujf\x) ^ g{x)U,{x)g^{x + ae,) 



(2.5) 



We use the following parametrization for the matrices U G SU(2): 
U = + iu ■ a = 



-U2 + iUi Uq — ms 



(2.6) 



where IL is the 2x2 identity matrix, the components of a = (cr^, o"^, cr^) are the Pauli 
matrices, «o G M, n G and ul + u.u = l. Therefore 



and from equations (|2.6| ) and (|2.7|) it follows that 

Tr [/ = Tr t/^ = 2uo . 
By using equations ( plTl) and ( |2.8|) we can also write 

U'^ = f/t = ^Tr f/ - f/ . 
\i V = vq + iv ■ a is another SU (2) matrix then, again using ( |2.6| ) and 



(2.7) 



(2.8) 



(2.9) 



we obtain 



-Tr ( U 



UqVq + U - V 
U ■ V , 



(2.10) 
(2.11) 



where the last step follows if we interpret a matrix U G SU{2) as a four-dimensional unit 
vector u = {uo,u). Finally, HUE SU{2), the matrix 



(2.12) 



* Here we do not consider the problem of searching for the absolute minimum of the minimizing function 
£ , which defines the so-called minimal Landau gauge [E8| . 
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belongs to the su (2) Lie algebra and is parametrized as 

A 



lu ■ a 



tUs U2 + lUi 

—U2 + iui —iu-i 



This matrix is traceless, anti-hermitian and 

Tr {AA^) = Ti{ -A^ 



2 u ■ u 



(2.13) 



(2.14) 



Let us now consider a one-parameter subgroup g{T; x) of SU{2) defined by 

g{T] x) = exp [ r7(x) ] 



(2.15) 



where the parameter r is a real number and the 7's are fixed elements of the su (2) Lie 
algebra given by 

7(3;) = i^{x)-a (2.16) 

with ^{x) E for all x. Then, for fixed {U^{x)}, the minimizing function £ can be 
regarded as a function of the parameter r, and its first derivative is given by the well-known 
expression 



£\0) 



dV 



E [^{x)]^ [(V ■ [x) 

X 

where the sum in the color index j is understood (j = 1, 2, 3) and 

[(V ■ A) {x)V ^-Y: [A,ix) - A,{x - ae,)] 



is the lattice divergence of 



[M^)], = ^Tr[A^(a;)a' 



(2.17) 



(2.18) 



(2.19) 



If {Ufj_{x)} is a stationary point of £{t) at r = {i.e. g{T,x) = lL,Va3) then we have 
^'(0) = for all {-fj{x)}. This implies 

[(V ■ A) (x)]^ = W X, J (2.20) 

which is the lattice formulation of the usual local Landau gauge-fixing condition in the 
continuum. By summing equation (|2.20|) over the components of x with /i 7^ z/ and using 



the periodicity of the lattice, it is easy to check that [0 if the Landau gauge-fixing condition 
is satisfied then the quantities 



X 
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(2.21) 



are constant, i.e. independent of Xy. From this it immediately follows that the longitudinal 
gluon propagator at zero three-momentum 



Duu{x.) ^ EEoTr(A.(a^)A.(0)) 
= ^Tr(g,(x,)v4,(0)) 



(2.22) 
(2.23) 
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is also constant. 

The numerical problem wc have to solve is therefore the following: for a given {i.e. fixed) 
thermalized lattice configuration {t/^(cc)}, we want to find a gauge transformation {g{x)} 
which brings the function 

= 1 - ^ E E Tr [ 9{x) U,{x) g\x + ae,) ] (2.24) 
^ ^ ^ n=i X 

to a local minimum, starting from a configuration g{x) = l_ for all x. In order to achieve 
this result it is sufficient to find an iterative process which, from one iteration step to the 
next, decreases the value of the minimizing function monotonically (descent methods). In 
fact, since S has a lower bound of (and an upper bound of 2), an algorithm of this kind is 

expected to converge. 

To find a simple iterative algorithm which minimizes S one may choose to update a single 
site variable g{y) at a time. In this case the minimizing function S becomes 

£[g(y)] = constant - ^p^Tr w(y) (2.25) 

where 

w{y) = 9{y)h{y) (2.26) 
and the "single-site effective magnetic field" h{y) is given by 

d 

h{y) = ^ [ UM 9\y + «e^) + Ul{y - ae^) g\y - ae^) ] . (2.27) 

The matrices h{y) and w{y) are proportional to SU (2) matrices, namely they can be written 
as 

h{y) = ^]deth{y) h{y) (2.28) 
w{y) = ^[d^^)w{y) (2.29) 

with h{y),w{y) e SU{2) and 

U{y) = ^Jdethiy) = ^Jdetw{y) . (2.30) 

Let us also define 

T{y) = TTw{y) = Ti w{y) / M{y) . (2.31) 

We want to consider the single site update g{y) — > g^'^'^'^\y) , which can alternatively 
be looked at as the multiplicative update 

9{y) - 9^''''^\v) ^ R^^^'^'^\y) g{y) , (2.32) 
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with _r("P'^"*'^) (y) G SU{2). Under this update, the variation of the minimizing function is 
given by 



AS{y) = _ Tr 



2dV 
2dV 



Tr 



9{y) - 9^''''"\y) 



Hy) 



(2.33) 
(2.34) 



To measure the length of the move g{y) — » 9^'^^'^\y) we can use| the quantity 

'9{y), 9^'''''"\y)] = V{y) ^ 



V 



i Tr \ [g{y) - g^^^^){y) ] [g{y) - g^^^^\y) ]^ 



2 - Tr [g^^''^\y) g\y)] 
2 — Tr _R("P'^"*^) (y) 



(2.35) 



which satisfies the defining properties of a distance function for any set of matrices and, if 
we interpret SU{2) matrices as four- dimensional unit vectors, it coincides with the standard 
euclidean distance in M'* [see formula ( p.ll| )]. 

In the next section the local quantities AS{y) and T^ly) will be used to illustrate 
the performance of the different local methods considered. In particular, their expressions 
will be written completely in terms of the tuning parameter for the algorithm (if needed), 
the square root Af{y) of the determinant of w{y) and h{y), and the trace T{y) of the 
normalized matrix w{y) . 



3 The Algorithms 

In this section we will describe the five algorithms for which we want to analyze the critical 
slowing-down. In particular, we will compare the implementation and performance of the 
four local algorithms^ considered in Sections p.l| - p.4| , by comparing their expressions for the 
quantities AS{y) and V^y) introduced in the previous section. 

As explained before, AS{y) measures by how much the "local energy" is reduced in a 
single step of the algorithm at site y, while T>{y) measures by how much the configuration 
at site y was effectively changed. Therefore they represent the two (possibly conflicting) 
tasks that a local algorithm is expected to perform: minimizing the energy at every site and 
at the same time moving efficiently through the configuration space. 

As we will see in Section |3.1| below, the Los Alamos method has the "best" possible value 
for AS{y) , i.e. it brings the "single site energy" to its absolute minimum in one iteration. 
We take the Los Alamos method as a basis for comparing all the other local algorithms we 



^ This choice is not useful in the case of global gauge transformations; in fact, for this kind of trans- 
formations we obtain a non-zero value for even though they do not really represent a move in the 
configuration space. 

^ A simple comparison of this kind is not possible for the Fourier acceleration method. 
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consider, which will typically have smaller A£{y) (in magnitude) and larger 'D{y) than the 
Los Alamos method, and will perform better. 

In order to make this comparison more quantitative we will also look at the ratios 



ne{y) = 



and 



(3.1) 



(3.2) 



for the various methods. They measure, respectively, the relative "loss" in minimizing £{y) 
and the relative "gain" in the length of the update, when compared with the Los Alamos 
method. In the next subsections these two quantities will be evaluated, for each local al- 
gorithm, as functions of M{y) and T{y) (defined in the previous section) and the tuning 
parameter. In particular, their limits as T{y) — > 2 will be computed. In fact, as discussed in 
Section |, these limits are useful to point out analogies between the algorithms and compare 
their efficiencies in fighting critical slowing-down. 



3.1 Los Alamos Method 

It is easy to see from equations (|2.33|) and ( p.34|) that the choice [1^, [13 

gineu^)^y) = g{LosAl.) ^y^ ^ ^y^ ^ 



or 



^{update) ^y^ _ ~t(^) 

gives the maximum negative variation of S: 



/\£^LosAl.)^y^ 



2dV 



:n?/)-2]<o 



(3.3) 
(3.4) 

(3.5) 



where 'T'{y) was defined in equation (p. 31 ). In other words, the move from g{y) to g^^"^^^-)[y'^ 
brings the function £ [giy)] to its unique absolute minimum. For this update we have, from 
equation ( |2.35| ), 

= J2 - T{y) . (3.6) 



3.2 Cornell Method 



Another possible choice for (^y) comes from considering an update of the form ( p.l5| ) 

with T = a and 7(?/) = —a^go (^V ■ A^^^^ (y). Then, from equation (|2.17|) we obtain 



[giy) 



-a 



0^ ' 
dV 



(3.7) 



and it is clear that the minimizing function decreases if a is small and positive. So we can 
define ^\ 

(3i 



j^iupdate)(^y^ = exp [-aa^go (v ■ A^^)^ (y) 
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Since we expect a? go (V ■ (y) to go to zero as the number of iterations increases, we 
can expand R^'^v'^"-^^) (^y^ to first order obtaining 



j^{update)^y^ oc [±-aa^ go (y ■ A^^^) (y) 



(3.9) 



here oc indicates that the matrix on the l.h.s. is proportional to the one on the r.h.s. 
(namely it has to be reunitarized) and the parameter a is a positive real number which has 
to be properly tuned, depending on the considered configuration, as discussed below. 
If we notice that the matrix w{y) [defined in (|2.26| )] satisfies the relation 



w{y)=w^{y) + 2a^go (V ■ A^^^) (y) 
we can rewrite equation (p.9|) as 



^(update) ^y-^ OC i + 



a 



w\y) 



w{y) 



and, by using equation]] ( ^l9|) with U = w\y) , we obtain 



j^{update) 



OC 



a 



y) 



IL + aw\y 



(3.10) 



(3.11) 



(3.12) 



Finally, reunitarizing (ty) and using ( p.32|) we have 



g{cornell)^y^ 



aN{y)w\y) + [1 - aN{y)r{y)/2] 1. 
l + a'M\y)[l-r\y)m 



9{y) 



(3.13) 



where w{y), M{y) and T{y) are defined respectively in ( p.29|) , (p.30|) and ( |2.31|) . In this 
case the variation of the minimizing function is given by 



^gicorneiD^y^ ^ <^'^{v) I 2«Ar(^) + [1 - aN{y) r{y)/2] T{y) 



2dV 



l + a^X^iy)[l-T\y)/A] 



(3.14) 



Since T{y) is in the interval [—2, 2] and a is positive this quantity is negative or zeroQ iff 
aM{y) G (0, 2]. Therefore the algorithm converges only if a is positive and small enough. 
On the other hand, if we evaluate the length of the move g{y) — >■ gi<^°'^"-^^'-) (^y^ -^e obtain 



from equation (2.35) 



-p{cornell)^y^ = ^2 1 - il + ft^^^ 



r\y) 



^1/2 \ 1/2 



(3.15) 



namely a should not be too small otherwise this length goes to zero. 



^ Note that equation (2.9) holds also for multiples of SU{2) matrices such as w^y) . 

^ To see this notice that A£ is zero at the end points T(y) = ±2 and negative for T{y) = 0, that there 
are no other zeros in this interval iff aN{y) € (0, 2] , and that for aN{y) > 2 the variation A£ approaches 
zero from above as T(y) goes to two. 
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It is easy to check that in the hmit of T{y) — 2 we get for the ratios TZs{y) and TZv{y) 
[defined respectively in and ( |3.2| )] 

^(corneu)^y^ [aAfiy) - if (3.16) 

and 

n^^orneiD^y^ - [a U {y) - 1] . (3.17) 

Therefore, as T{y) approaches 2 we have that, at least for the case 1 < aN'{y) < 2, 
the gain in the length of the move with respect to the Los Alamos method is linear in 
[aN'{y) — 1] , while the loss in the minimizing of the energy is quadratically small. This 
illustrates why the algorithm will perform better than the Los Alamos method. For further 
discussion, see Section 5. 



3.3 Over relaxation Method 



The standard overrelaxation method ||T6[ is a local algorithm in which, instead of using the 
update 

-^^'^^''''■\y) = w\y) g{y) (3.18) 



described in Section 3.1, we use the substitution 



9{y) 



w\y) 



9{y) 



(3.19) 



where the overrelaxation parameter uo varies in the interval 1 < u; < 2 and has an optimal 
value which is volume- and problem-dependent. Of course, for a; = 1, we have g^°''^'\y) = 



9 



(LosAl.) , 



y) while, for = 2, we obtain 



^(--)(^) = w^{y)w\y)g{y) 
and it is easy to check that [see equation ( |2.34| )] 



2dV 



Tr 



w[y) 



w\y) 







(3.20) 



(3.21) 



namely for uo = 2 the algorithm does not converge, as the energy is never decreased. 
Finally, for 1 < a; < 2, we can write 



w\y) 



(3.22) 



Therefore we can interpret this update as a move from g{y) to 

g{over)(^y^ "passing through" 

g{LosAi.)^yy |.j^jg way, the minimizing function £ [g{y)] will not go to its absolute minimum, 
even though its variation AS will still be negative. For computing w\y) one uses the 
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binomial expansion]^ 



w\y) 



T(uj + 1] 



oo 

„5 n\r{u + l 



n] 



1. 



(3.24) 



Since the matrix w''{y) is expected to converge to the identity matrix IL, this series can be 
truncated after a few terms, followed by reunitarization of the resulting matrix; for example, 
if only two terms are kept, we have 



± + U 



■9[y) 



'l + uj{LU-l)[2-T{y) 
In this way the variation of the minimizing function is given by [ see equation ( |2.34| ) 

{l-uj)T{y) + 2uj 



(3.25) 



AS^-\y) = { ny) 



l + uj{uj-l)[2-T{y) 



(3.26) 



Since T{y) G [—2, 2] and uo G (1, 2) it is easy to check that this variation is always negative 
or zeroQ. The length of the move [formula ( |2.35| )] is given by 



V^°'"^''\y) 



2-uj[2- T{y) 



1/2 



l + u{uj-l)[2-T{y) 



(3.27) 



As an illustration of the improvement with respect to the Los Alamos method, let us 
consider oj slightly larger than 1 . Expanding the expressions for A£(°^^^) and around 
a; = 1, we obtain 



A8^-'\y) = my) - 2] { 1 - 1 (.; - 1)^ [2 + r{y) ] + O {{u - if) } 



and 



V^-'\y) = ^2-T{y) { 1 + ^ (c. - 1) [2 + T{y) ] + O {{u - if) } 



^ We could also write the matrix w{y) as 

w{y) = 1 cos 7 + in ■ a sin 7 = exp (i"fn ■ a) 



(3.28) 



(3.29) 



(3.23) 



with 7 G [~7r, tt) , rt G and n ■ n = 1. Then [ (y) would be given by exp (—ijLun- a) where the 
product 7u; should be considered modulo 2 7r so that 70; G [— 7r,7r) . In any case we are interested in the 
limit in which w'^{y) approaches the identity matrix JL, name ly 7 — > . By expanding exp(— i7wn- a) 
around 7 = 0, and reunitarizing we obtain again formula ( 3.25 ). 



Furthermore, it can be proved that, if a; > 1/2, the variation A£ is zero only at T{y) = 2 while, if 
u! < 1/2 , this happens at both end points T(y) = ±2 . Then it is easy to check that A£ is always negative 
or zero if cj G [0, 2] . 
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Figure 1: Plot of the ratios TZf^^ (y) and 71^^'^^ (y) as functions of Tiw = T , for 
comparison between the overrelaxation method at a; = 1.9 and the Los Alamos method. 

namely the correction with respect to the variation ( |3.5| ) is positive and quadratic in (u; — 1), 
while the correction with respect to the length ( |3.6| ) is positive and linear in {u — 1). Thus, 
already for a value of uj slightly larger than one, what we lose in minimizing S, compared 
with the Los Alamos method, is "smaller" than what we gain in the length Vly) of the move 
for the update and therefore the relaxation process should be speeded up. 

More generally, these features can be seen from the behavior of the quantities TZ^^^'^^^ (y) 
and 7i}^^^\y) [defined respectively in ( |3.1| ) and ( |3.2| )]. As an example, we plot in Figure 
1 these two ratios as a function of T{y) and with uj = 1.9 . In particular, in the limit 
'^{y) ^ 2 we obtain 

n'-r'^\y) {uj ~lf (3.30) 

and 

n^-'\y) - (^-1) . (3.31) 

It is interesting to note that the behavior is qualitatively the same as the one for the Cornell 
method, discussed at the end of the previous subsection. 



3.4 Stochastic Overrelaxation Method 



The stochastic overrelaxation method is also a local algorithm. In this case, instead of 
always applying the descent step g{y) g(^°^^'--)(^y'^ one uses the new update 



9iy) 



g^''°'\y) 



w\y) 



g{y) with probability p 



(3.32) 



9 



(LosAl.) 



{y) 



with probability 1 — p 



with < p < 1. Of course for p = this algorithm coincides with the Los Alamos method 
while, for p = 1, it does not converge at all since, as we saw in formula ( p.21| ), the value of 
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the minimizing function £ [g{y)\ remains constant. However, for p G (0, 1), the fact that, 
with probabihty p, a big move is done in the configuration space without changing the value 
of £ [g{y)] has, again, the capability of speeding up the relaxation process. To check this 
point we can compute the length of the move 



g{y) w\y) 



9[y) 



(3.33) 



by using equation 



we can rewrite 



w^'iy) 



as 



w\y) 



w 



[y) T{y) 



1. 



and [see equation (|2.35|) ] we easily obtain]^ 



2 + T{y) p(^°^^'-)(z/) . 



(3.34) 



(3.35) 



Roughly speaking, we can say that the stochastic overrelaxation method alternates up- 
dates that give the maximum negative variation of £ with steps that produce "very long" 
moves in the configuration space, without increasing the value of the minimizing function. 
This is similar in spirit to the idea behind the so-called hybrid version of overrelaxed algo- 
rithms 
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which are used to speed up Monte Carlo simulations with spin models, 
lattice gauge theory, etc. In these algorithms, n microcanonical (or energy conserving) up- 
date sweeps are done followed by one standard local ergodic update (like Metropolis or 
heat-bath) sweeping over the lattice. Actually, for the gaussian model, it has been proven 
p5| that the best result is obtained when the n microcanonical steps are picked at random, 
namely when n is the average number of microcanonical sweeps between two subsequent 
ergodic updates. This is essentially what is done in the stochastic overrelaxation method, 
with n/ {n + 1) equals in average to p or, equivalently, n equals in average to p/ {1 — p). 
Finally, formula ( |3.34| ) can also be used in order to reduce the overhead of the algorithm; 



we can in fact rewrite ( ^.32| ) as 



(stoc) , 



y) 



h^{y)T{y) — g{y) with probability p 
h'^{y) with probability 1 — p 



(3.36) 



where h was introduced in ( p.28 ) ], namely instead of computing matrix products of the 



form h'^{y) g\y) h^{y) we just have to do a simple linear combination of h^{y) and g{y) 



3.5 Fourier Acceleration 



The idea of Fourier acceleration fl^ is very simple. If we consider the Cornell method, it 
is immediate from formula (|3.8|) that its convergence is controlled by the quantity g^ (V ■ 



It is also straightforward to check that for this algorithm, as T(y) goes to two, the ratios TZsiy) and 
T^viy) , defined at the beginning of Section 3, are both equal to one, with probability p, and to zero, with 
probability 1 — p. 
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A^^^){y) . For the abelian case in the continuum it can be shown |T^, by analyzing the 



relaxation of the different components of this matrix in momentum space, that 

(V ■ A^3^)t (fc) ^ (V ■ A^3))o (k) exp (-a p\k) t ) , (3.37) 

namely each component decays as exp{—ap'^{k) t ) , where t indicates the number of sweeps. 
This means that their decay rates are approximatively equal to l/(ap^(fc) ) . Therefore, if we 
choose a oc p^ax obtain that the slowest mode has a relaxation timef^ r proportional to 
Pmax/Pmin ■ ^ lattice with points on each side we have p^ax ^(1) oc (9(A^~^) , 

namely 

r oc A^^ _ (^3 38) 

From this analysis it is clear how, for the abelian case, the relaxation process can be speeded 
up: given the matrix gQ^V ■ A^^^){x) , we take its Fourier transform, we multiply each 
component in momentum space by Pmax/P^i^) > evaluate the inverse Fourier transform 
and, finally, the result is used in equation ( |3.9|) instead of the original matrix qq (V ■ 
a more concise form we can write: 



j^iupdate) (^) 0, ^ _ ^-1 J ]^ F[aa^go {V- A^^)^] \ (y) ; (3.39) 



p2(fc) 

where F indicates the Fourier transform and is its inverse. In this way we should obtain 
that the components in momentum space of go (V ■ A^^^){y) decay as 

exp ^) = exp {-apl^j) (3.40) 

which, with the choice a oc p"^^ , gives r oc 0{1) for every component. 

Of course, for the non-abelian case, this analysis becomes more complicated. Neverthe- 
less, it is still believed[3 that the Cornell method will have ||, |14| the behavior (|3.38|) , and 
that the strategy to be used in the Fourier acceleration is given by the modified update ( |3.39| ). 



The main difficulty arises from the fact (see again |14]) that, instead of the eigenvalues of 
the laplacian {i.e. instead of p^(fc) ), we have to consider the eigenvalues of the operator 
d ■ D, where D is the covariant derivative. Thus, the relaxation time r will be proportional 
to the ratio of the largest over the smallest eigenvalue of d ■ D and the eigenvectors of this 
operator should be used to decompose the divergence of A^^\ This is of course not easy to 
be implemented in a numerical simulation and, therefore, the eigenvectors of the laplacian 
are used also in the non-abelian case. The hope is that the non-abelian nature of the fields 
does not make the behavior of d ■ D , in momentum space, too different from that of the 
laplacian. Actually this is more than just hope since, in the lattice Landau gauge, the link 



variables {U^{x)} are fixed as close as possible to the identity matrix (see ||28|, Appendix 



For an exact definition see formula (4.11). 

Note that the behavior ( 3.38 ) corresponds to dynamic critical exponent 2 for the Cornell method. This 
is in contradiction with the analysis given in Section 5, which predicts an exponent z « 1, by analogy with 
the overrelaxation method. Our results (see Section 7) corroborate the latter prediction. 
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A) and therefore the operator d ■ D should be, in some sense, a "small modification" of the 
laplacian. 



The practical implementation of the Fourier acceleration is also quite simple. In fact, 
we have to evaluate go (V ■ A^^^){x) at each lattice site and then use formula (|3.39|) — 
where now F has to be interpreted as a standard Fast Fourier Transform subroutine [|l^ — 
to find Ri^P'^'^t^) (y) at the Efiven lattice site. Of course, to reduce the number of times the 
FFT is used, a checkerboard update should be employed. For our FFT subroutine we used 
as a basis in momentum space the functions exp(2 7rifc ■ x) , where k has components 
given by 

afc^iV^ = 0,1, ...,iV^-l (3.41) 

and fi = 1, . . . ,d . In this case the eigenvalues of the laplacian operator are given by the 
well-known formula 

= 4 E sin2(7raA;J (3.42) 



and the largest eigenvalue p^ax obtained when 



L 2 ^ 



(3.43) 



for all fi = 1, . . . ,d . 

Finally, it is important to observe that formula (|3.39| ) is singular when p^{k) is zero. 
However, the zero-frequency mode of the divergence of A^^"* does not contribute to the update 



( pIOD ; in fact, by using the periodicity of the lattice and formula (|2.18|) , it is easy to check 
that 

a2^oE(V-^(^^)(^) = 0. (3.44) 



X 



Thus, in equation (|3.39| ) when p^(fc) is equal to zero, we can set the value of Pmax/p'^i^) 
any finite value without affecting the performance of the method. 



4 Critical Slowing-Down 



To check the convergence of the gauge fixing 0, |I^ |T6|, |18|, ^ several quantities have been 
introduced: 



V 



X j=l 



1 



ei(t) = £{t-l)-S{t) 
62 (t) 

64 (i) 



EE [(v-A)(^) 



— 5] - Tr ] 1 - R(^P<^^t^) (x) 1_ - R(^pd'^'^) (x) 
^ X 



max 

X 



1 - -Tri?("P'^'^*")(a;) 



1 _ ^ E Tri?("P'^'^*^)(aj) 



(4.1) 
(4.2) 

(4.3) 

(4.4) 

(4.5) 



X 
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where t indicates the number of lattice sweeps and, when not indicated, the expressions on 
the r.h.s are evaluated after t sweeps of the lattice are completed. All these quantities are 
expected to converge to zero exponentially and with the same rate [§] even though their sizes 
can differ considerably. Actually, it is easy to see that 



2a° 



E 

X 



Tr J^("P'^''*'^) (x) 



2 65 ^ 2 64 . 



(4.6) 



By using equation (|3.10| ) we can also rewrite 62 as 



62 



V 



X 



1 - 



THx) 



(4.7) 



We know that, if the algorithm converges, the matrix /^("P'^"*^) (aj) should approach the iden- 
tity matrix IL as the number of iteration increases or, equivalently, that its trace should be 
very close to two at large t . This implies [ see the expressions of Ri'^P'^'^^^) (^x) for the various 
local methods: formulae (P^), ( |3.13| ), ( |3.25| ) and ( |3.32| ) ] that w'^{x) IL , and therefore 
T{x) is also very close to two0. Therefore if 64 is of order e ^ 1 then also 63 , 65 and 62 
should be of this order; taking into account these relations, we decided to look only at the 
quantities ei, 62 and 64 . 

We also expect that the spatial fluctuations of the quantities Qu{x^), defined in (|2.21| ), 
and of the longitudinal gluon propagators, defined in ( p.22|) , go to zero exponentially. Indeed, 
the fact that Qu{x^) should be constant is being increasingly used as a check of the accuracy 
of the gauge fixing [p3l To check this more precisely, we introduced a new quantity 



defined as 



U=l 



3K 



where 



3 

■E E 

j=l x„=l 
1 



[a 



(4i 



(4.9) 



For each of these quantities, in the limit of large t , we can introduce a relaxation 
time Ti by the relation 

ei{t) ^ a exp{-t/Ti) , (4.10) 



namely 



Ti 



lim — I 



'ei{t + l)/ei{t)] ■ 

As said before, we expect all these relaxation times to coincide and be equal to r. 



(4.11) 



14 



From equations ( ^.34 ) and ( ^.35 ) it is clear that also A£{x) and V^x) go to zero. In particular, since 64 
is proportional to maxx [^{x)]^ , it is obvious that "D^x) goes to zero exponentially. This tells us that, when 
the condition T{x) < 2 is satisfied (usually after a few sweeps) , the algorithm "moves" very slowly through 
the configuration space and therefore improving the accuracy of the gauge fixing becomes very costly. 



17 



To analyze the critical-slowing down of an algorithm we have to measure r for different 
pairs^ of N and f3 = A/g^, but at "constant physics". This means that we have to keep 
the ratio N/C, constant, where the correlation length ^ is given by the inverse of the square 
root of the string tension k, i.e. ^ = The string tension for two-dimensional SU{2) 

lattice gauge theory (in the spin-| representation) is given, in the infinite volume limit, by 

log^€. (4.12) 



hi(3) 



where is the modified Bessel function. Thus we have 



/i(/3) 



(4.13) 



which, in the limit of large /3 , gives 




A/3 



+ 0{(3- 



(4.14) 



Therefore a constant ratio N/S, is equivalent, in this limit, to keeping the ratio N"^ / (3 fixed. 
The values for the pairs (A^, (3) have been chosen in such a way that N ^ ^] thus the finite 
size effects should be negligible. All the pairs {N,j3) used are reported^ in Table |l|. In the 
same table we report the value of the corresponding correlation length ^ obtained by using 
equation (|]T|). We have chosen N'^/p = 32 and N/^ ^ 7. 



Once all these values of r are obtained, we can try a fit of the form 



T 



(4.15) 



and find the dynamic critical exponent z for that algorithm. It is important to recall that the 
value of z obtained in this way is independent of the "constant physics" chosen [N"^ / (3 = 32 
in our case). On the contrary, this is not the case for the constant c. In particular we expect 



the relaxation time r, and therefore c, to increase as the link couplings U^{x) in ( ^.24 ) 
become more "random", 
ratio N'^ / (3 increases). 



z.e as {3 decreases for a given lattice size (and the value of the 



5 Tuning of the Algorithms 

The implementation of all the algorithms considered in this work — except for the Los 
Alamos method — requires the tuning of a parameter: a for the Cornell method and the 
Fourier acceleration, uj for the overrelaxation method and p for the stochastic overrelaxation. 

From now on we always consider lattices with N^^ = N for all /i = 1, . . . , c? . 

For the case of the Fourier acceleration we considered N — 8, 16,32,64. Note that we restricted our 
lattice sizes to powers of 2, because of the way in which the Fast Fourier Transform subroutine we used 
is designed. (The application of this routine is not limited to these lattice sizes, and it can be modified to 
work with general values of N , but the use of powers of 2 makes it most efficient.) 



18 



N 


8 


12 


16 


20 


24 


28 


32 


36 


64 




2.0 


4.5 


8.0 


12.5 


18.0 


24.5 


32.0 


40.5 


128.0 




1.09 


1.65 


2.24 


2.83 


3.42 


4.00 


4.58 


5.16 


9.22 



Table 1: The pairs {N,j3) used for the simulations and the correlation length ^ evaluated 
in the infinite volume limit. 



algorithm 




Ky) 


Los Alamos 


1 





Cornell 




aAf{y)T{y)/2 - 1 


overrelaxation 




LO- 1 


stochastic overr. 


T{y) with probability p 
1 with probability 1 — p 


1 with probabihty p 
with probability 1 — p 



Table 2: The coefficients a{y) and h{y) for the four local algorithms considered in this 
paper. 



This is, of course, a potential disadvantage of all these methods and makes the study of 
their critical slowing-down more difficult: in fact, for each pair {N, j3) , the value of the 
parameters should be tuned in such a way that the value of r is minimized. This is usually 
done heuristically since, as explained in the Introduction, no rigorous analyses are available 
for these algorithms. However, analytic estimates for the optimal choice of u are indeed 



known for the overrelaxation method applied to other numerical problems 0, [IT], |2J, ^ 
in all these cases, in the limit of large lattice size A^, it has been found that 

= 1 + L,t/N ^^-'^ 

where the constant Copt is problem dependent. This result is usually adopted as a guess 
|l6| , p!9| also for the optimal choice of uj when the overrelaxation method is applied to Landau 
gauge fixing. 

In order to obtain simple formulae like (|5.1| ) for the optimal choice of a (in the case of 
the Cornell method) and p, we decided to compare the four local algorithms considered in 
this work. It is interesting to notice that they can all be defined by the update|^ 

9{y) - g'-''''"^ oc a{y)h^{y) - b{y) g{y) (5.2) 

where the coefficients a{y) and b{y) are given in Table |^. Moreover, we see from our simu- 
lations that in all methods, usually after a few sweeps, we have T{y) < 2 . This is evident 
by looking at the decay of 62 (t) (see Figures |^ and ^ for typical examples) and considering 
formula ([4.7|). Actually, since the gauge-fixing procedure is stopped when 62 (t) is smaller 



This is, of course, not surprising if we notice that equation (5.2) is the most general linear local update 
we can introduce to minimize the function A£'-'""'''\y). 
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than 10 (see end of Section |^), the condition T < 2 is satisfyed for the most part of our 
simulations. With this in mind we can write 

g{y) - ^7^"^'") ^ ~a{y)h\y) - Ky)g{y) (5.3) 

where the coefficients a{y) and h{y) are obtained from the coefficients a{y) and h{y) by 
imposing the condition T{y) = 2. 

This simple analysis seems to suggest (see Table 2) that the Cornell method is equivalent 
to the overrelaxation method if we make the substitution 



a 



M{y) u . (5.4) 



The same substitution is suggested by considering the ratios TZsiy) and 7lx>{y) defined in 
( ^.11 ) and (|3.2|); in fact, as T{y) goes to 2, we have that the formulae for these ratios for 
the Cornell method and for the overrelaxation method coincide if the above substitution is 



employed (see end of Sections |3]^ and ^]3| respectively). 

If this analysis is correct we should obtain for the Cornell method the same dynamic 
critical exponent as for the overrelaxation method, i.e. z 1 . This will be verified in 
Section 7. As a further check of this "equivalence" between the two methods we can compare 
the optimal choices for their tuning parameters, obtained in our simulations. We notice that 
while uj and a are fixed parameters throughout the run, N'{y) changes with the iterations, 
and we are interested in its value as T{y) — > 2 . Moreover, N{y) is a local quantity, and 
therefore we consider its space average, which can be easily estimated in this limit. In order 
to do this, let us rewrite the minimizing function as 

£ = l-T^EETr (f/^^H^) + K^n^-ae,)]' 

= l-J^ET. ^E^WrW (5.5) 

which, in the limit T{y) 2, gives 

£-1-5^I:A^W- (5.6) 

In other words, the space average of Af{x) is given by 

2dil-£stat) , (5.7) 

where Sgtat is the value of the minimizing function at the stationary point. Of course this 
value is not known a priori but, for fixed values of (3 and l^, its order of magnitude can be 
easily estimated with just a few numerical tests. Using this result, we are able to make a 
numerical comparison between the tuning parameters for the two methods (see Section 7), 
finding very good agreement. 

One may also attempt to establish a relation between the overrelaxation and the stochas- 
tic overrelaxation. For example, we can try to write the update ( p.32|) as an average of the 



two cases with weights p and 1 — p obtaining (in the limit T{y) 2) 

g^''"'^\y) ^ {l+p)h\y) -pg{y), (5.8) 
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which suggests the substitution 

p ^ {to - I) . (5.9) 



However, if we now look at the ratios Tigiy) and IZviy) for the update ( p. 81 ) we obtain, as 
T{y) goes to 2, 

nf'"'\y) - p (5.10) 

and 

<°^^(?/) -p; (5.11) 



the second formula, if compared to ( p. 311) , seems to be consistent with the substitution (|5^ ) 
while the first [compared to (|3.3CI|) ] suggests the relation 

p ^ {uj - if . (5.12) 

These two possibilities will be tested numerically (see Section ^ by plotting the quantities 
{uj — 1) / p and {uj — 1 Y / p as a function of A^. 

Finally, we do not hazard here any hypothesis on the tuning of a for the Fourier accel- 
eration method. 



6 Numerical Simulations 

To thermalize the gauge configuration {U^{x)} at a fixed value of the coupling /3 , we used 
a standard heat-bath algorithm [^. In order to optimize the efficiency of the code, we 
used two different SU (2) generators (methods 1 and 2 described in Appendix A in , with 
hcutoff = 2). 

With the pairs (A^, /3) that we considered (see Section we should always have ^ ^ 
and therefore we expect all the temporal correlations to decay exponentially. As a check we 
measured, for all the pairs (iV, /?), the integrated autocorrelation tim^ for the Wilson loop 
iy(l, 1) and for the Polyakov loop P (indicated respectively as Tint,Wi '^mt,p)- More- 
over, for the pairs (A^, f3) used for the Fourier acceleration method, we also measured the in- 
tegrated autocorrelation time Tint^Wi for the Wilson loops W{l,l) with I = 2,4,8, ...,N/2. 

In practice, we started all our runs with a random {U^{x)} configuration and we did 
5000 sweeps for thermalization. The configurations used for gauge fixing were separated 
by 100 sweeps, in order to get a statistically independent sample. After discarding 4900 
sweeps out of a total of 54900 sweeps, we evaluated Tint,Wi 'Tint,p by using a window of 
4 Tint , which is a reasonable choice if the decay is roughly exponential. For the observables 
we considered we obtainedj^ 0.5 ^ Tint ~ 10 . Noticing that Tint = 0.5 indicates that the 
data are uncorrelated, we can conclude, as expected, that the system decorrelates rather 

See 1^ for a definition of integrated autocorrelation time and for a description of tlie automatic 
windowing procedure used to measure it. 

From our data it is clear that, for a fixed lattice size, the Wilson loop with size I « ^ has the largest 
integrated autocorrelation time among the quantities we considered. Nevertheless, we obtain Tint i$ 2.5 
for all our lattice sizes, except for 64^ (used only for Fourier acceleration), where we get Tint,Wi ^ 10 ^'^^ 
I = 8,16. 
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fast and that the configurations used for testing the gauge-fixing algorithms are essentially 
statistically independent. 



The tuning of the parameters a;, p and a — respectively for the overrelaxation, the 
stochastic overrelaxation, the Cornell and the Fourier acceleration methods — was done very 
carefully. More precisely, we divided it in three parts. In the first step, we considered a few 
values of the parameter spread in a large interval. For example we used u = 1.1, 1.2, . . . , 1.9 
or p = 0.1, 0.2, . . . , 0.9 . For each of these values we gauge fixed 10 different configurations, 
measured all the 's and averaged the results. In this way we were able to select a smaller 
interval for the parameter (usually of length ^ 0.2 for the overrelaxation or the stochastic 
overrelaxation methods) which was used in the second step of the tuning. In this case 20 
configurations were analyzed for each value of the parameter (and these values were typically 
separated by 0.01 for the overrelaxation and stochastic overrelaxation methods). In the last 
step, the length of the interval was further reduced and 100 configurations were gauge-fixed 
for each value of the parameter. 

A total of 201.7 hours of CPU were used for the four methods requiring tuning. Of these, 
21.6 were used for the first level of tuning, 36.1 for the second level and the remaining 144 
for the third level. 

For the Los Alamos method no tuning is needed, but we found that, in this case, the 
fiuctuations for the relaxation times are larger and therefore more configurations (500 , to 
be exact) had to be considered, for a total of 23.3 hours of CPU. 

Finally, for measuring the relaxation time (with i = 1,2,4 and 6) we did a chi- 
squared fitting of the functions logej(t) which, if relation (|4.10|) is satisfied, should be a 



straight line. Indeed this was usually the case already after a few initial sweeps of the 
lattice, at least for the quantities ei , 62 and 64 . On the contrary, the decay of eg is really 
smooth and monotonic only for the Fourier acceleration method. As an example we show, 
in Figure H, the behavior of 62 and eg as functions of t for the Cornell and the Fourier 
acceleration method. We also show, in Figure ^, the decays of the four quantities for the 
stochastic overrelaxation method. 
We use the condition 

62 < 10-^2 (6.1) 

to stop the gauge fixing, in order to ensure that enough data are produced for the fitting and 
that, when the procedure is stopped, essentially only the slowest mode has survived. To get 
rid of the initial fiuctuations, we used only the last 100 data when the total number of sweeps 
Nsw was larger than 200 , or the second half of the data when less than 200 sweeps were 
necessary to fix the gauge. For z = 1 we have also to take into account possible fiuctuations 
of ei{t) around zero, which appear when the minimizing function is fixed within the machine 
precision. Therefore, for this quantity, we also discarded the last 50 sweeps, if Ng^ > 200, 
or the last one quarter of the data if A^^^ was smaller than 200. 

7 Results and Conclusions 

Our final data for the relaxation times are reported, for the different methods, in Tables 
We show for each algorithm only the relaxation time T2 for the quantity 62 defined in 
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Cornell Method 




50 100 150 

t 



Fourier Acceleration 




10 20 30 40 50 60 

t 



Figure 2: Plot of the decays of the quantities and eg as functions of t for (a) the Cornell 
method, with a — 0.481, and (b) the Fourier acceleration method, with a. — 0.160. Both 
plots were done for \& lattices. 
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50 100 

t 



150 



Figure 3: Plot of the decays of the quantities ei , 62 , 64 and eg as functions of t for 
the stochastic overrelaxation method, with p = 0.75 on a 16^ lattice. The two almost 
superposed curves are 62 (t) and 64 (i) {e2{t) is the "smoother" curve). 
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Indeed, we checked that, for all methods and pairs (iV, /?), the four measured relaxation 
times were in agreement within error bars. We also show the optimal choice for the tuning 
parameter (when needed), the number of sweeps necessary to reach the stopping condition 
( |6.1D , and, for the Cornell method, the value of the minimizing function (used in Section 
7.3 for comparison with the overrelaxation method). Averages are taken over the different 
configurations that were gauge fixed for each pair (A^, f3) . 

7.1 Critical Exponents 

We now proceed to the evaluation of the dynamic critical exponents z . In Tables |To|-[l^ we 
present the results of the fits to the ansatz ( [4.15| ) for T2 for the various methods. We do a 
weighted least-squares fit in several "steps" , discarding at each step the values of smaller 
than Nmin ■ In this way we try to rule out some of the smaller values of A^ as finite-size 
corrections; this is very important since we are dealing with very small lattice sizes. We do 
so for all possible values of Nmm , and decide which one gives the best fit for z by comparing 
and confidence levels for the different Nmin 's. As can be seen from our tables, these 



finite-size corrections are negligible already at lattice sizes 12 or 16. In Figure ^ we plot 
together the values of T2 and the fitting curve for our preferred fit for the various algorithms. 

Our results for the dynamic critical exponents are in agreement with what is generally 
expected, namely we find: z pa 2 for the Los Alamos method, z ^ 1 for the overrelaxation 
|19| and the stochastic overrelaxation method, and z ^ for the Fourier acceleration 



method. For the Cornell method, as mentioned in Section p.5| , a simple analysis (based 



on the abelian case in the continuum) would give z ^ 2 , and this is what is generally 



believed P, On the other hand, our comparative analysis between the Cornell and the 



overrelaxation methods in Section 5 would suggest z ^ 1 . As can be seen from Table 11 
our results show the latter behavior. Furthermore, in Section 7.3 we will verify the relation 
between the tuning parameters for these two methods found in Section 5. 

Actually, the dynamic critical exponent z for the Cornell method is even slightly smaller 
than one. The good performance of this method could be understood by noticing that the 
value oiaN'lx) changes during the gauge-fixing process. In particular, with a few numerical 
tests, it can be checked that the space average value of Af{x) increases with the iterations. 
Moreover, as we will see below, the final value of aAf{x) is equal to the optimal value found 
for uj in the overrelaxation method. So, in a sense, we have an overrelaxation method whose 
parameter u ^ aAf{x) increases with the iterations^ from an initial value Uq up to the 
asymptotic value Uopt ■ It is well known [|l^, 22 1 that in overrelaxed algorithms the optimal 



strategy is precisely to vary the parameter u from an initial value 1 to a larger asymptotic 
value Uopt , which is usually done by using Chebyshev polynomials. It is conceivable that the 
Cornell method does this variation "automatically" and this could explain why it performs 
slightly better than the overrelaxation method. 



It should be stressed that this scenario is very quaUtative. In particular, as explained in Section |5|, the 
relation between cu and aM{x) is established only when T < 2, i.e. it should not be used for the initial 
sweeps of the lattice. 
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filgorithni 


fi 




Los Alamos 


0.2445 ± 0.0008 


0.5197 ±0.0113 


Cornell 


0.1490 ± 0.0054 


4.5212 ±0.5858 


overrelaxation 


0.1061 ±0.0106 


3.2357 ±0.3482 


stochastic overr. 


0.0151 ±0.0039 


0.9642 ±0.0214 


Fourier acceleration 


0.0483 ±0.0148 


1.0219 ±0.0374 



Table 3: The ratios ri and r4 for each algorithm. Averages are taken first over the gauge- 
fixed configurations and then over the different pairs {N, (3) . The error bars are one standard 
deviation. 

7.2 Checking the gauge fixing 

For each gauge-fixed configuration we also measured the ratios 



with 2 = 1,4 and 6 . For the cases i = 1 and i = A this quantity is essentially independent 
of the configuration and of the lattice size; therefore, for each algorithm, after averaging over 
all the configurations, we take a final average over all the pairs (A^, (3) . The results are given 
in Table ^. From that it is clear that the quantities ei , 62 and 64 not only decay with the 
same rate (as said above) but also have the same order of magnitude. 

The situation for the ratio ee is considerably more complicated. In fact, its value depends 
strongly not only on the algorithm and on the lattice size N but also on the underlying 
configuration^ {f/^(a;)}. Namely, this quantity fluctuates so much that if the average is taken 
over all the gauge-fixed configurations, at a fixed lattice size, the corresponding standard 
deviation is often comparable in magnitude to the average itself. As an example, see Table 
^ where we show our results for all the methods on a 32^ lattice. From these data (see also 
Figures ^ and ^ it is clear that the Fourier acceleration method achieves a much faster decay 
for eg than the Los Alamos method, the Cornell method and the overrelaxation method. 
Actually this was expected. In fact, by using one of these three local methods it can be 
easily checked that, even when the condition (|6.1|) is satisfied, the quantities Qy [defined 
in (|2.21|) ] are usually not constant but show a kind of long-wavelength spatial fluctuation^. 
The Fourier acceleration method treats all the wavelength in the same way and therefore it is 
not surprising that it is very effective in reducing these spatial fluctuations. More surprising 
is the good performance of the stochastic relaxation method which, although a local method, 
also appears to be very efficient in reducing the fluctuations of Q,j . Why this happens is 
not clear to us. 



^-^ That this should be the case was, somehow, expected since ee is a less "local" quantity th an ei , 
and 64 , and therefore it represents a more sensible check for the lattice Landau gauge condition ( 2.20| ). 

See also Figure 12 in O and Figure 1 in 
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8-lgoritlini 


^6 


inin tq 


nicLX tq 


Los Alamos 


10.7 X 10«±4.9 X 10" 


4939.7 


2.23 X 10^ 


Cornell 


19.3 X 10*^ ± 19.0 X 10« 


86.3 


1.90 X 10*^ 


overrelaxation 


2.4 X 10^ ± 1.2 X 10^ 


34.1 


1.1 X 10^ 


stochastic overr. 


3.9 X 10^ ± 2.0 X 10^ 


4.93 


1.8 X 10^ 


Fourier acceleration 


5.6 X 10^ ± 4.5 X 10^ 


0.84 


4.5 X lO'' 



Table 4: The ratio rg (average, minimum and maximum value) for each algorithm. Averages 
are taken over the gauge-fixed configurations for = 32 (for other lattice sizes the results 
are similar). The error bars are one standard deviation. 

7.3 Discussion of the Tuning of the Algorithms 

Let us now discuss the tuning of the different methods. The values for the optimal choice of 
the various parameters, for different pairs {N,(3) , are reported in Tables |^-|^. An estimate 
of their uncertainties is also indicated. From our simulations we noticed that a good tuning 
becomes more and more important as the lattice size increases. At small lattice sizes, in fact, 
the relaxation time r displays a kind of plateau around the minimum while, as A^ increases, 
the absolute minimum becomes more and more pronounced. The uncertainties indicated in 
these tables are therefore slightly under-estimated for the smaller lattice sizes, and slightly 
over-estimated for the larger lattice sizes. In Figure ^, as an example, we show typical graphs 
of our tuning parameters for some of the larger values of A^ , done at the "third level" of 
tuning (see Section 6). 

We now try to verify, using our data, the expressions suggested in Section 5 for the 
various tuning parameters. 



In order to check the relation (|5 . 1| ) for the optimal choice of the overrelaxation parameter 
u , we rewrote that equation as 

(7.2) 



'^opt _ Copt 



uJopt N 

and fitted our results to find a value for the constant Copt ■ After discarding the datum for 
A^ = 8, we obtained Copt = 1-53 ± 0.35 . In Figure |^ we show, together, the data and the 
fitting curve. 

For the parameter a of the Cornell method we do not have a simple formula as ( |7.2| ) but 
we have conjectured, in Section H, a relation between ujopt and Uopt- Namely we suggested 
[see formulae (^.4|) and ( p.7|) ] 

^^opt = O^opt 2d {1 - Sstat ) , (7.3) 

where d is the dimension of the lattice and Sgtat is the value of the minimizing function at 
the minimum. By using the data reported in Tables |^ and |^ we plotted together, in Figure 
^, both sides of this equation. The agreement is clearly good. 

Finally we checked the two relations introduced in Section ^ between uj and the tuning 
parameter p for the stochastic overrelaxation method. In particular we plotted, in Figure |^ 
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the two ratios 

Popt Popt 

as a function of . If one of these relations [ see formulae ( ^.91) and (|5.12|) ] is correct 



we should obtain, for the corresponding ratio, a constant value 1. From our data it is not 
possible to reach a definite conclusion, but the first hypothesis, namely 

^^^^^^ = 1 , (7.5) 

Popt 

seems to be slightly better verified. 

7.4 Computational Cost of the Algorithms 

To check the computational cost of the algorithms we estimated the CPU time Tgf necessary 
to update a single site variable g{x) by using the fortran function mclock. As expected, 
the four local methods have very similar values for Tgf and essentially independent of the 
volume. In particular we found Tgf ^ 9/is for the Los Alamos method, Tgf ^ 9.5/is for 
the Cornell method, Tgf ^ 11.5/is for the overrelaxation method and Tgf ^ 10.5/is for the 
stochastic overrelaxation method. 



For the case of the Fourier acceleration method Tgf should increase [|T0| as log A^ . We 
did a least-squares fit of our data to the ansatz 

Tgf = A\ogN + B (7.6) 

and we obtained A = 48.91 ± 0.11 and B = —54.42 ± 0.45, both measured in microsec- 
onds. In Figure 1^ we show the points and the fitting curve. For our range of lattice sizes, Tgf 
for the Fourier acceleration varied from 56 us to 148 us . Of course, this "loss" in efficiency, 
with respect to the local algorithms, has to be taken into account when the computational 
cost of this algorithm is analyzed. In fact, even though the Fourier acceleration method 
succeeds in eliminating critical slowing-down, and therefore is more advantageous than the 
improved local method in terms of the number of sweeps A^^^ required to achieve gauge 
fixin^^, its performance is effectively better only at very large lattices. 

To illustrate this, let us compare the true CPU time needed to gauge fix a configuration 
using Fourier acceleration and our best improved local algorithm: the Cornell method. Since 
we want to evaluate the total CPU time we have to look at the number of sweeps A^^^^, , needed 
in average to achieve gauge fixing, as a function of A^. This quantity behaves in a manner 
similar to the relaxation time r , namely it diverges with A^ as a power of some dynamic 
critical exponent Znsw ■ This exponent should be similar to, but strictly smaller than the 
exponent z for the relaxation times. In fact, A^s^ is a quantity involving the behavior of 
the whole gauge-fixing process, and therefore it includes faster modes (modes with "smaller 
exponents") for the first few iterations, while the relaxation time r is evaluated in the 
limit of large t, i.e. when essentially only the slowest mode has survived. The exponents 



^•^ From Tables ||-^ it is clear that the number of sweeps Nsw increases roughly linearly with the lattice size 
for the improved local methods, while it remains essentially constant for the Fourier acceleration method. 
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Overrelaxation Method, N = 20 

8 . I \ \ \ \ \ \ \ \ \ \ \ 



7.5 



7.0 



6.5 



6.0 
1.80 



1.85 



1.90 



Stochastic Overrelaxation, N 

16 



15 



28 



14 



13 



12 



11 



10 
0.80 



0.82 0.84 



0.86 
P 



01 



0.90 0.92 



Figure 4: Plot of the third level of tuning respectively for (a) the overrelaxation method, at 
lattice size 20, and (b) the stochastic overrelaxation method, at lattice size 28. Error bars are 
one standard deviation. Note that points are correlated, since the same 100 configurations 
are used for each value of the tuning parameter. 
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0.02 



0,01 I — III < < < < — < — II 

5 10 20 50 100 

N 



Figure 5: Plot of the ratio (2 — uJopt)/oJopt (symbol: o) as a function of A^. The solid line 
is the fitting curve Copt/N with Copt = 1-53. 



Relation between oj and a 



10 20 30 40 



Figure 6: Plot of ujopt (symbol: o) and UoptM (symbol: □) as a function of with M 
given by equation (|5.7| ). 
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Relation between cj and p 




Figure 7: Plot of (cUopt — l)/Popt (symbol: o) and (cUopt — l)^/Popt (symbol: □) as a function 
of N. In order to check the two hypotheses introduced in Section 5, the constant curve 1 is 
also shown. 



Znsw for the various methods can be obtained, together with the respective proportionality 
constants Cnsw , from a fitting of the data in Tables ||-^ analogously to what was done for 
the z's. For the Cornell method we obtain Znsw ~ 0.77 (which should be compared to 
z ~ 0.82 for the relaxation time) and Cnsw ~ 18 , while for the Fourier acceleration we get 
Znsw ~ (as expected) and Cnsw ~ 67. From these estimates, the value Tgf ^ 9.5/is for 
the Cornell method and the fit (|7.6| ) for the Fourier acceleration method, we can get the 
following approximate expressions for the time, in microseconds, necessary to gauge fix a 
configuration: 

r 9.5 (18A^°-^^) N^fis Cornell method 

CPUtime ^ I (7.7) 
[ ( 49 log — 54 ) 67 A^^ jjs Fourier acceleration 

In Figure ^ we show a plot of these two functions (divided by the volume A^^). Clearly, 
Fourier acceleration becomes the method of choice only at lattice sizes of order of 350 !! 
Of course this analysis is very machine- and code- dependent. In particular we remark that, 
at the present stage, our code for the Cornell method has been considerably optimized, while 
the one for the Fourier acceleration can still be improved, hopefully increasing its running 
speed by a constant factor, which we think could be as high as 2. Moreover, if we used a 
condition on the quantity eg to stop the gauge fixing, instead of the one given in equation 
(p.l|), then the computational cost of the Cornell method would increase much more than 



that of the Fourier acceleration method, as is clear from the discussion in Section \l.2[ In 
any case, it seems unlikely that the Fourier acceleration method would become the method 
of choice at lattice sizes smaller than around 100 sites. 
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Fourier acceleration; N^^^ = 8 

200 I — , , , ^ — r- 




' — III < < < < — < — II 

5 10 20 50 100 

N 



Figure 8: Plot of CPU time Tgj as a function of N for the Fourier acceleration method. 
The solid line represents a least-squares fit to the ansatz Tgf = A log N + B , with A = 
48.91 ± 0.11 and B = -54.42 ± 0.45 (both in microseconds). 
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Figure 9: Comparison between CPU times (divided by volume) as a function of the lat- 
tice size for the Cornell and the Fourier acceleration methods. The "almost" straight fine 
corresponds to the Cornell method. 
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7.5 Conclusions 

From our numerical simulations it is clear that the Fourier acceleration method is very 
effective in reducing critical slowing-down for the problem of SU{2) Landau gauge fixing in 
two dimensions. On the other hand, its computational cost is much larger than that of the 
improved local methods, and therefore an accurate analysis should be always done to decide 
which method to use. The result of this analysis, as stressed in the previous subsection, will 
depend on the code, on the machine and on the condition used to stop the gauge fixing. 
From our data it seems that, at least up to lattice size of order 100, the improved local 
methods should be always preferred. 

From the point of view of computational cost, the Cornell method is clearly the best 
among the improved local methods. However, if the condition used to stop the gauge fixing 
is not ( |6 . 1|) , this conclusion could be different. In particular we saw that the stochastic 
overrelaxation is very effective in relaxing the value of the quantity cq , defined in ( [4.8[ ). 

All the improved methods, including the Fourier acceleration, have the disadvantage of 
requiring tuning. However, the relations for the overrelaxation, the Cornell and the stochastic 
overrelaxation methods that were introduced in Section 5 and checked in Section 7.3 make 
the tuning of these methods simpler. Moreover, as reported in Section 6, the values of r 
for the improved methods [for a fixed pair {N,f5)] are much more stable than for the Los 
Alamos method. Therefore, in order to find the optimal choice of tuning parameter within 
a few per cent, it suffices to perform a few numerical tests. 

Finally, from the discussion in Sections 4 and 7.2, it is clear that the quantities ei-es 
are essentially equivalent as a check of the goodness of the gauge fixing. Namely, when one 
of these quantities is measured, the evaluation of any of the others does not provide any 
new information. On the contrary, the quantity eg represents a more sensible check of the 
gauge-fixing condition and, in our opinion, should always be evaluated. 

In a future paper we will try to extend this analysis to the more interesting case of 
lattice gauge theory in four dimensions. 
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N 


r 


sweeps 


8 


14.71 ±0.68 


330 ± 16 


12 


29.87 ± 1.14 


600 ± 18 


16 


53.32 ±2.00 


1054 ± 34 


20 


86.55 ±3.17 


1650 ±91 


24 


125.01 ±7.32 


2159 ±79 


28 


157.22 ±5.56 


2797 ± 103 


32 


205.21 ±8.03 


3372 ± 105 


36 


282.78 ± 14.29 


4389 ± 155 



Table 5: The relaxation time and the number of sweeps for the Los Alamos method. Error 
bars are one standard deviation. 



TV 


T 


a 


sweeps 


niin. fimc. 


8 


3.22 ±0.30 


0.489 


86 ±5 


0.1171 ±0.0012 


12 


4.54 ±0.33 


0.484 


117±6 


0.0688 ± 0.0005 


16 


6.31 ±0.63 


0.481 


152 ±9 


0.0430 ±0.0003 


20 


6.92 ± 0.40 


0.480 


181 ±8 


0.0296 ± 0.0002 


24 


8.60 ± 0.53 


0.482 


217 ± 8 


0.0218 ± 0.0001 


28 


9.52 ± 0.61 


0.483 


232 ±9 


0.0168 ±0.0001 


32 


10.70 ±0.43 


0.484 


260 ±7 


0.0134 ±0.0001 


36 


11.65 ±0.44 


0.484 


278 ±7 


0.0108 ±0.0001 



Table 6: The relaxation time, the coefficient a , the number of sweeps and the value of the 
minimizing function for the Cornell method. Error bars are one standard deviation. The 
estimated error on the parameter a is about 0.002 for all lattice sizes. 



N 


r 


UJ 


sweeps 


8 


2.77 ±0.35 


1.72 


77 ±9 


12 


3.56 ±0.15 


1.77 


95 ±3 


16 


4.84 ±0.20 


1.83 


136 ±5 


20 


6.51 ±0.32 


1.86 


208 ± 32 


24 


7.84 ±0.38 


1.88 


208 ±9 


28 


9.04 ±0.52 


1.90 


249 ± 11 


32 


10.34 ±0.62 


1.90 


257 ±9 


36 


12.33 ±0.61 


1.92 


306 ±7 



Table 7: The relaxation time, the coefficient uj and the number of sweeps for the overrelax- 
ation method. Error bars are one standard deviation. The estimated error on the parameter 
UJ is about 0.01 for all lattice sizes. 
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N 


T 


P 


sweeps 


8 


3.35 ±0.28 


0.64 


96 ±5 


12 


4.38 ±0.18 


0.72 


129 ±4 


16 


6.48 ±0.46 


0.78 


189 ± 11 


20 


7.84 ±0.30 


0.84 


234 ±7 


24 


9.47 ±0.46 


0.85 


281 ± 13 


28 


11.28 ±0.67 


0.88 


331 ± 13 


32 


12.67 ±0.53 


0.89 


366 ± 12 


36 


14.95 ± ().<S7 


0.90 


421 ± 15 



Table 8: The relaxation time, the coefficient p and the number of sweeps for the stochastic 
overrelaxation method. Error bars are one standard deviation. The estimated error on the 
parameter p is about 0.02 for all lattice sizes. 



N 


T 


a 


sweeps 


8 


3.08 ±0.33 


0.193 


80 ±7 


16 


3.30 ±0.34 


0.170 


84 ±7 


32 


3.11 ±0.21 


0.160 


92 ±8 


64 


3. 40 ± 0.33 


O.KiO 


93 ± 9 



Table 9: The relaxation time, the coefficient a and the number of sweeps for the Fourier 
acceleration method. Error bars are one standard deviation. The estimated error on the 
parameter a is about 0.003 for all lattice sizes. 
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N 

^ * rmn 


z 


c 


X 


o 




9441 + 09"^^ 




12 


1.986 ± 0.042 


0.2174 ± 0.0284 


4.443 (5 DF, level = 48.751 %) 


16 


1.965 ± 0.060 


0.2332 ± 0.0448 


4.196 ( 4 DF, level = 38.017 %) 


20 


1.919 ± 0.090 


0.2727 ± 0.0810 


3.718 ( 3 DF, level = 29.353 %) 


24 


2.030 ± 0.171 


0.1857 ± 0.1082 


3.131 ( 2 DF, level = 20.900 %) 


28 


2.281 ± 0.238 


0.0779 ± 0.0636 


0.819 ( 1 DF, level = 36.551 %) 


32 


2.722 ± 0.543 


0.0164 ± 0.0313 


0.000 ( DF, level = 100.000 %) 



Table 10: Weighted least-squares fit for r = cN^ at A^^/ j3 = 32, using lattice sizes > Nmm, 
for the Los Alamos method. Errors represent one standard deviation, and "DF" stands for 
"degrees of freedom" . Confidence level is the probability that would equal or exceed the 
observed value. The line in boldface marks our preferred fit. 



N ■ 


z 


c 


x' 


8 


0.854 ± 0.049 


0.5509 ± 0.0877 


1.134 (6 DF, level = 98.002 %) 


12 


0.849 ± 0.062 


0.5611 ± 0.1156 


1.114 (5 DF, level = 95.284 %) 


16 


0.825 ± 0.089 


0.6088 ± 0.1833 


0.977 (4 DF, level = 91.329 %) 


20 


0.859 ± 0.106 


0.5409 ± 0.1940 


0.608 (3 DF, level = 89.456 %) 


24 


0.762 ± 0.168 


0.7605 ± 0.4401 


0.045 (2 DF, level = 97.759 %) 


28 


0.788 ± 0.283 


0.6936 ± 0.6878 


0.032 (1 DF, level = 85.749 %) 


32 


0.721 ± 0.471 


0.8809 ± 1.4627 


0.000 (0 DF, level = 100.000 %) 



Table 11: Weighted least-squares fit for r = c at N'^ / j3 = 32, using lattice sizes N > Nmin, 
for the Cornell method. Errors represent one standard deviation, and "DF" stands for 
"degrees of freedom" . Confidence level is the probability that would equal or exceed the 
observed value. The line in boldface marks our preferred fit. 



^ min 


z 


c 


x' 


8 


1.094 ± 0.045 


0.2390 ± 0.0329 


3.381 ( 6 DF, level = 75.974 %) 


12 


1.120 ± 0.048 


0.2205 ± 0.0325 


1.062 ( 5 DF, level = 95.745 %) 


16 


1.120 ± 0.069 


0.2202 ± 0.0485 


1.061 (4 DF, level = 90.034 %) 


20 


1.063 ± 0.108 


0.2673 ± 0.0949 


0.577 ( 3 DF, level = 90.157 %) 


24 


1.101 ± 0.164 


0.2340 ± 0.1296 


0.479 ( 2 DF, level = 78.696 %) 


28 


1.239 ± 0.301 


0.1445 ± 0.1513 


0.185 ( 1 DF, level = 66.707 %) 


32 


1.491 ± 0.659 


0.0590 ± 0.1375 


0.000 ( DF, level = 100.000 %) 



Table 12: Weighted least-squares fit ioT t — c a,t N'^ / P — 32, using lattice sizes > Nmin, 
for the overrelaxation method. Errors represent one standard deviation, and "DF" stands 
for "degrees of freedom" . Confidence level is the probability that would equal or exceed 
the observed value. The line in boldface marks our preferred fit. 
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^ min 


z 


c 




8 


1.048 ± 0.043 


0.3395 ± 0.0449 


3.673 ( 6 DF, level = 72.089 %) 


12 


1.086 ± 0.050 


0.3002 ± 0.0465 


1.330 (5 DF, level = 93.178 %) 


16 


1.034 ± 0.081 


0.3566 ± 0.0940 


0.680 ( 4 DF, level = 95.377 %) 


20 


1.065 ± 0.098 


0.3217 ± 0.1029 


0.355 ( 3 DF, level = 94.929 %) 


24 


1.084 ± 0.171 


0.3017 ± 0.1751 


0.338 ( 2 DF, level = 84.455 %) 


28 


1.115 ± 0.329 


0.2702 ± 0.3084 


0.325 ( 1 DF, level = 56.842 %) 


32 


1.407 ± 0.609 


0.0966 ± 0.2062 


0.000 ( DF, level = 100.000 %) 



Table 13: Weighted least-squares fit for r = cN"" at N"^/ (3 = 32, using lattice sizes > Nmin, 
for the stochastic overrelaxation method. Errors represent one standard deviation, and "DF" 
stands for "degrees of freedom" . Confidence level is the probability that would equal or 
exceed the observed value. The line in boldface marks our preferred fit. 



Nrnu, 




c 


A" 


8 


0.036 ± 0.064 


2.8513 ± 0.6062 


0.751 (2 DF, level = 68.703 %) 


16 


0.040 ± 0.102 


2.8080 ± 1.0094 


0.748 ( 1 DF, level = 38.712 %) 


32 


0.157 ± 0.169 


1.8020 ± 1.1286 


0.000 ( DF, level = 100.000 %) 



Table 14: Weighted least-squares fit for r = c at N'^/ /3 = 32, using lattice sizes N > Nmin, 
for the Fourier acceleration method. Errors represent one standard deviation, and "DF" 
stands for "degrees of freedom" . Confidence level is the probability that would equal or 
exceed the observed value. The hne in boldface marks our preferred fit. 
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FOURIER ACCELERATION; N 




Figure 10: Power-law fit for (a) tlie Los Alamos method, (b) the Cornell method, (c) 
the overrelaxation method, (d) the stochastic overrelaxation method and (e) the Fourier 
acceleration method. 
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